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A Unitary Fermi Supersolid: The Larkin-Ovchinnikov Phase 
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We present strong theoretical evidence that a Larkin-Ovchinnikov (loff/fflo) pairing phase is 
favoured over the homogeneous superfluid and normal phases in three-dimensional unitary Fermi 
systems. Using a Density Functional Theory (dft) based on the latest quantum Monte-Carlo 
calculations and experimental results, we show that this phase is competitive over a large region 
of the phase diagram. The oscillations in the number densities and pairing field have a substantial 
amplitude, and a period some 3 to 10 times the average interparticle separation. Within the dft, the 
transition to a normal polarized Fermi liquid at large polarizations is smooth, while the transition 
to a fully-paired superfluid is abrupt. 

PACS numbers: 03.75.Ss, 03.75.Hh, 05.30.Fk, 



The BCS mechanism for fcrmionic superfluidity is 
rooted in the notion of pairing: Can superfluidity survive 
in polarized systems with unequal numbers? This ques- 
tion remains largely unanswered, even though it is funda- 
mental to many forms of matter, including superconduc- 
tors, nuclear matter, and high density QCD. The prospect 
of observing exotic polarized supcrfluids has been revived 
in two-component cold-atoms gases with s-wave interac- 
tions, especially in the unitary limit where the scattering 
length diverges: \a\ — > c». Here, the physics is universal, 
and inherently strongly coupled, depending solely on the 
densities. These systems are experimentally tenable (see 
Ref. [1] for a review), and exhibit a remarkable diversity 
of polarized phases. 

Clogston and Chandrasekhar [2] noted that the normal 
phase competes with BCS superfluidity when the chemical 
potential difference between the species becomes compa- 
rable to the energy gap. Kohn and Luttinger [3], how- 
ever, showed that interactions render Fermi surfaces un- 
stable at sufficiently low temperatures, suggesting pair- 
ing of higher angular momenta. This effect is exponen- 
tially suppressed in weak coupling, but may be strong 
enough in unitary gases to support symbiotic p-wave 
superfluids [4]. Another proposal by Fulde and Fer- 
rell (ff) [5], and Larkin and Ovchinnikov (lo) [6] - 
anisotropic/inhomogeneous polarized superfluids, widely 
referred to as loff or fflo states - have also been vig- 
orously sought (see [7] for reviews), but firm results have 
been sparse: Experimentally there have been claims of 
quasi-two-dimensional fflo states in heavy-fcrmion su- 
perconductors [8], but no 3d realizations have been re- 
ported. Other proposals include deformed Fermi sur- 
faces [9] and gapless (breached pair) superfiuids [10]. 

We present here strong evidence that an inhomoge- 
neous Larkin-Ovchinnikov (lo) state [6] may be realized 
in cold polarized unitary Fermi gases. Our approach is 
novel in several respects: 1) it is the first calculation to 
find a completely self-consistent LO solution in three di- 
mensions; 2) the calculation is based on a Density Func- 
tional Theory (dft) incorporating the best Monte-Carlo 



calculations and measurements of the unitary Fermi gas; 
and 3) includes both pairing and self-energy correlations. 

Previous calculations of LO states have not been fully 
self-consistent, often relying on approximate forms of 
spatial variations, or uncontrolled Ginzburg-Landau ex- 
pansions (see e.g. Refs. [7, 11, 12]). Furthermore, self- 
consistent treatments are typically based on mean-field 
or Bogoliubov-de-Gennes (BdG) calculations, which do 
not properly account for many-body effects such as the 
Gorkov-Melik-Barkhudarov corrections [13] that lead to 
significant decreases in the pairing gap. Finally, most cal- 
culations account for only the pairing condensation en- 
ergy, which is exponentially suppressed in weak-coupling, 
while the LO state has density variations that can signif- 
icantly change the unsuppressed normal correlation en- 
ergy ( "Hartree" terms) . Mean-field and BdG calculations 
neglect these crucial correlation contributions: Without 
them, LO states are not competitive at unitarity. 

According to the theorems of Hohenberg and Kohn, 
a dft exists for any system of fermions. At unitarity, 
the structure of the functional is strongly constrained 
by dimensional arguments, and thus its determination is 
greatly simplified. The remarkable accuracy of this ap- 
proach for symmetric systems - as demonstrated in Ref. 
[14] - gives us the confidence to extend the approach to 
polarized systems. To model the polarized unitary Fermi 
gas, we use an asymmetric (aslda) generalization of the 
superfiuid local density approximation (slda) employed 
in Ref. [14], expressed in terms of the following densities 
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where M„(r), w„(r), and En are the quasiparticle wave- 
functions and energies [14]. We use the same functional 
form as Ref. [14] , but allow the parameters to depend on 



the local asymmetry x{r) = ni,{r)/na{r). The resulting 
ASLDA energy density £{r) has the form 
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where aa{r) = a[a;(r)] and ai,{r) = Q![l/.T(r)] are the 
inverse effective masses in units of m^^ defined in terms 
of the single function a(x), (3(x) = (3(1/ x) parametrizes 
the normal interaction, and {ua + nt)^'^ /j = l/^cff + A 
defines the effective coupling ^eff that is regulated with 
the cutoff A as described in Refs. [14, 15]. The forms of 
a{x) and f3{x) are well constrained by Monte-Carlo data, 
as described below and in Ref. [15]. The equations for 
the quasiparticle wave- functions u„(r) and Vn{r) follow 
by minimizing the grand-canonical functional 
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where iiafi are the chemical potentials corresponding to 
the two fermion species, and V{r) is the local pressure. 

We must now specify the forms of a{x) and I3(x). First 
we analyze the symmetric superfluid phase as described 
in Refs. [14, 15]. By matching the Monte-Carlo val- 
ues [16, 17, 18, 19] for the parameters ^ = £sf/£fg — 
0.40(1), 77 = A/ef = 0.504(24), and the single quasipar- 
ticle excitation spectrum, we determine the inverse effec- 
tive mass a(l) « 1.09(2), the constant 7-1 = -0.091(8), 
and the implied energy of the symmetric normal phase 
Cat = £n/£fg = a(l) + /?(!) = 0.57(2). In Ref. [14] 
it was erroneously stated that q;(1) could be extracted 
from the values of ^ and 77 alone. A more careful analysis 
shows that the quasiparticle dispersion [17] must also be 
fit, resulting in the modified SLDA parameters presented 
here. The inverse masses aa,b are also known for the 
fully-polarized gas. The majority species is unaffected, 
whereas the minority species has the effective "polaron" 
mass m/m* « 1.04(3) [20, 21], constraining the end- 
points a{0) = 1 and a(oo) = 0.96(3). To determine 
the function f3{x), we consider the energy of the inter- 
acting normal state, setting A = 0. This has been well 
constrained by Monte-Carlo calculations [16, 20], and, 
along with the parametrization of a{x), a fit to this data 
uniquely specifies the function (3{x) (see Fig. 1). 

In this paper, we consider only the simplest LO states, 
with spatial modulations in a single direction (z). Unlike 
the FF state [5], the lo state [6] does not break time- 
reversal invariance, and thus belongs to a different sym- 
metry class, as was emphasized by Yoshida and Yip [11]. 
We do not consider ff states here as they are typically 
not competitive with the LO states, (ff states break 
time-reversal invariance, and require additional terms 
in (2) to restore Galilean invariance [14, 15].) 
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FIG. 1: (Color online) The dimensionless convex func- 
tion g{x) [22] that defines the energy density £{na,ni,) = 

11^(6^^)^''^ [""5(^)]^''^- The points with error-bars (blue 
online) are the Monte-Carlo data from Refs. [20]. The fully- 
paired solution g{l) = (2^)'^'^ is indicated to the bottom 
right, and the recent MIT data [23] is shown (light x) for com- 
parison. The phase separation discussed in Refs. [20] is shown 
by the Maxwell construction (thin black dashed line) of the 
f^'-order transition y = yN-sF in Fig. 2. The LO state (thick 
red curve) has lower energy than all pure states and phase 
separations previously discussed. The Maxwell construction 
of the weakly l°*-order transition y = yLo-sp in Fig. 2 is 
shown by the thick dashed line (red). 



The self-consistency equations are solved by discretiz- 
ing the Hamiltonian along z with a discrete variable rep- 
resentation (dvr) basis [24] of period L, while integrating 
over the perpendicular momenta and the Bloch states. 
Our quasiparticle wave- functions thus satisfy the condi- 
tions Un{x,y,z + L) = e*'^'"u„(a;,2/,z), Vn{x,y,z + L) = 
e"^"u„(x,?/, z), and are plane waves in the {x,y) plane. 
We minimize the truncation error due to the finite DVR 
basis set by using a smoothed version of the hybrid strat- 
egy [25] , summing discrete states with energies less than 
a cutoff Ec, while integrating over the remaining higher- 
energy semi-classical states. 

We start by specifying chemical potentials fXa and /ib, 
and an ansatz for A(z) ex sin(27rz/L) containing a node 
at z = 0, and then use a Broyden iteration scheme [26] 
to find a self-consistent solution. The choice of basis and 
iteration preserve the node at z = 0, converging to either 
a self-consistent LO state, or degenerating to a homoge- 
neous normal state with A(z) ~ everywhere. 

The resulting self-consistent states depend on the ex- 
ternal parameter L. To find the physical LO state, we 
vary L to find the spontaneously chosen length scale 
L = Llo that minimizes the potential f2 (maximizes the 
average pressure V). The search is aided by the relation- 
ship LdV/dL = 2£ — ZV between L, the average pres- 
sure, and the energy density [15], ensuring the unitary 
relationship T' = |5 is satisfied by the physical state. 

At unitarity, one may fully characterize all stable 
phases by the single parameter y = ^ih/ Ha as described 
in Ref. [22]. We use the grand canonical ensemble, where 
only pure phases appear, and which properly accounts 
for the phase separation that occurs at T**-order transi- 
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FIG. 3: (Color online) A single LO period showing the spatial 
dependence of the pairing field A(z) (top) and the number 
densities of the majority (dotted) and minority (solid) species 
(bottom) at the values of y G {yLO-N,yLO-SF) marked by x 
in Fig. 2. Units are fixed in terms of /i_ [28]. 



FIG. 2: (Color online) The dimensionless convex function 
h{y) [221 that defines the average pressure Vdj-a, jJ-t) ~ 
Ht^J [Ma/i(y)l^''V(67r^) is constrained to the thin dotted 
triangular region [22]. The interacting normal state pres- 
sure [20] defining our ASLDA functional constrains this fur- 
ther (thin blue line), and displays a f^'-order transition at 
yiv-sF where normal and SF phases could coexists (Maxwell 
construction in Fig. 1). The LO state has an even higher pres- 
sure (thick red line), replacing much of this region, including 
the former yN-sF transition. The y dependence of the ampli- 
tude of the pairing field A — max{|A(2:)|} and the period L 
are shown inset. Sample profiles for the states marked x are 
shown in Fig. 3. Units are fixed in terms of /i_ = (/la ~ P^b)l'i- 
[28]. 



tions (kinks in h{y), discontinuities in h'{y): see Fig. 2). 
We start by describing the homogeneous and isotropic 
states supported in the ASLDA functional: For y < yo [22], 
the system is a fully-polarized non-interacting Fermi gas 
{Na); between yo < y < yN-sF the highest pressure cor- 
responds to a partially polarized two-component Fermi 
gas; and above j/jv-sf < y < 1, the fully-paired superfluid 
(sf) has the highest pressure. The point y^-sF, where 
the pressures of the partially polarized normal and fully- 
paired superfluid states are equal, is where the phase- 
separation discussed in Rcfs. [20] would occur. Here, the 
competition to LO from the normal and superfluid states 
is minimized, and the LO state is most likely to occur. 
For 2/ > 1, the picture is reversed with the species a «-> 6 
exchanged. Our ASLDA parametrization does not ad- 
mit any stable homogeneous gapless superfluid (breached 
pair) states [10]. 

As shown in Fig. 2, we find competitive LO solutions 
for a large range of the parameter y e {yLo-N,yLo-sF) 
with finite periods in the range Llq-n ^ L < Llq-sf- 
At j/to-N, the transition appears to be second order, with 



max{|A(z)|} — > vanishing smoothly from the LO phase 
to the normal phase, while at yLo-sF, the order param- 
eter abruptly looses its spatial oscillations at a finite pe- 
riod Llo-sf- Due to the presence of a node in A(z), the 
only possibility for a smooth transition here would be for 
the period to diverge Llo-sf ^^ oo, thus this transition 
appears to be weakly first order. The remaining normal 
states between yo < y < y^o-N would be susceptible to 
the Kohn-Luttinger instability, and are candidates for the 
symbiotic p-wave supcrfluids discussed in Refs. [4]. To 
study this possibility requires an extension of the ASLDA. 

Fig. 3 shows the typical structure of a LO state. The 
pairing amplitude increases smoothly from zero at j/io-w, 
where the profile is almost sinusoidal, to a critical value 
slightly less than Aq at y^o-sF, while the minority com- 
ponent exhibits large oscillations that break translation 
invariancc, giving the LO state the crystalline properties 
of a quantum solid. The majority component exhibits 
much smaller oscillations because the larger local kinetic 
energy density suppresses gradients. These fiuctuations 
induce large oscillation in the mean-field potentials (not 
shown), and have a significant impact on the normal cor- 
relation energy. For this reason, all the terms in the en- 
ergy density functional are critical for a proper descrip- 
tion of the LO phase. 

If no other phases compete, one should observe that the 
LO-SF transition coincides with the termination point yi 
of the partially polarized phases PPa [22]: yLo-sF = J/i- 
A more complicated crystalline LO state with modula- 
tions in all three directions may further increase the av- 
erage pressure, making yi > yLo-sF- Current errors of 
the current Monte-Carlo calculations and experiments 
do not allow us to distinguish between these two cases. 
These results are summarized in Fig. 2, where it can be 
seen that these periodic LO solutions occupy a substantial 



portion of the phase diagram, and lead to a significant 
increase in the average pressure. The Legendrc trans- 
formed results have been included in Fig. 1 to facilitate 
comparisons with the Monte-Carlo data [20] and recent 
experiments [23]. 

In conclusion, we have shown that the ASLDA pro- 
vides a valuable tool for quantitatively evaluating in- 
homogeneous phases. By incorporating the latest non- 
perturbative data about unitary Fermi gases, we have 
presented strong evidence that a new form of matter, 
such as a crystalline LO phase, is waiting to be found 
in the partially polarized regime of cold unitary Fermi 
gases. This would be the first example of a Fermi super- 
solid at unitarity, and with the large pairing gap, there 
is a good chance of successfully studying this state with 
cold atoms. In the experiments to date, the shells where 
LO phases may exist are too thin to allow for a complete 
LO period L. However, traps can be adjusted so that 
the LO phase will occupy a larger spatial region, allow- 
ing for several loff oscillations to occur. Unlike LO in 
weak coupling, the amplitude of the density fluctuations 
in the minority component is large and comparable to 
that of vortices [27] in unitary gases. This will provide 
the most direct signature of unitary supersolid matter, 
and an clean way to study the loff phase. 
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